d <-#read_csv(paste0(home, "/data/for-analysis/dat.csv")) %>% read_csv("https://raw.githubusercontent.com/maxlindmark/perch-growth/master/data/for-analysis/dat.csv") %>%filter(age_ring =="Y", # use only length-at-age by filtering on age_ring!area %in%c("VN", "TH"))
Rows: 338460 Columns: 12
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (6): age_ring, area, gear, ID, sex, keep
dbl (6): length_mm, age_bc, age_catch, catch_year, cohort, final_length
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
filter: removed 43,886 rows (13%), 294,574 rows remaining
# Minimum number of observations per area and cohortd$area_cohort <-as.factor(paste(d$area, d$cohort))d <- d %>%group_by(area_cohort) %>%filter(n() >100)
# Minimum number of observations per area, cohort, and aged$area_cohort_age <-as.factor(paste(d$area, d$cohort, d$age_bc))d <- d %>%group_by(area_cohort_age) %>%filter(n() >10)
# We can also consider removing individuals where the SE of k is larger than the fitIVBG %>%#mutate(keep = ifelse(k > quantile(k_se, probs = 0.95), "N", "Y")) %>%mutate(keep =ifelse(k < k_se, "N", "Y")) %>%#filter(row_id < 10000) %>%ggplot(aes(row_id, k, ymin = k_lwr, ymax = k_upr, color = keep)) +geom_point(alpha =0.5) +facet_wrap(~keep, ncol =1) +geom_errorbar(alpha =0.5) +NULL
mutate: new variable 'keep' (character) with 2 unique values and 0% NA
# Add back cohort and area variablesIVBG <- IVBG %>%left_join(d %>%select(ID, area_cohort)) %>%separate(area_cohort, into =c("area", "cohort"), remove =TRUE, sep =" ") %>%mutate(cohort =as.numeric(cohort))
Joining with `by = join_by(ID)`
left_join: added 2 columns (area_cohort_age, area_cohort)
> rows only in x 0
> rows only in y (146,393)
> matched rows 142,023 (includes duplicates)
> =========
> rows total 142,023
mutate: converted 'cohort' from character to double (0 new NA)
# Compare how the means and quantiles differ depending on this filteringIVBG_filter <- IVBG %>%drop_na(k_se) %>%#filter(k_se < quantile(k_se, probs = 0.95)) %>% filter(k_se < k)
group_by: 2 grouping variables (cohort, area)
summarize: now 306 rows and 6 columns, one group variable remaining (cohort)
ungroup: no grouping variables
# No difference at all really. We should probably just discuss that with this model, achieving biologically reasonable parameter values and a good fit to data are sometimes two different things. In our case, we just want a representative value of the growth (as in time to reach average maximum size in the population) that accounts for the entire growth trajectory of an individual, for each area and cohort.
Rows: 380 Columns: 4
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (2): area, model
dbl (2): year, temp
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
rename: renamed one variable (cohort)
VBG <- VBG %>%left_join(pred_temp, by =c("area", "cohort"))
left_join: added 2 columns (temp, model)
> rows only in x 0
> rows only in y ( 74)
> matched rows 306
> =====
> rows total 306
# Save data for map-plotcohort_sample_size <- IVBG %>%group_by(area, cohort) %>%summarise(n =n()) # individuals, not samples!
group_by: 2 grouping variables (area, cohort)
summarise: now 306 rows and 3 columns, one group variable remaining (area)
VBG <-left_join(VBG, cohort_sample_size, by =c("cohort", "area"))
left_join: added one column (n)
> rows only in x 0
> rows only in y ( 0)
> matched rows 306
> =====
> rows total 306
write_csv(VBG, paste0(home, "/output/vbg.csv"))# Calculate the plotting order (also used for map plot)order <- VBG %>%ungroup() %>%group_by(area) %>%summarise(min_temp =min(temp)) %>%arrange(desc(min_temp))
ungroup: no grouping variables
group_by: one grouping variable (area)
summarise: now 10 rows and 2 columns, ungrouped
order
# A tibble: 10 × 2
area min_temp
<chr> <dbl>
1 SI_HA 7.83
2 BS 6.22
3 BT 5.87
4 FM 5.87
5 SI_EK 5.49
6 MU 5.16
7 FB 5.04
8 HO 3.99
9 JM 3.37
10 RA 3.06
write_csv(order, paste0(home, "/output/ranked_temps.csv"))nareas <-length(unique(order$area)) +2# to skip the brightest colors that are hard to readcolors <-colorRampPalette(brewer.pal(name ="RdYlBu", n =10))(nareas)[-c(6,7)]
Plot VBGE fits
# Sample 30 IDs per area and plot their data and VBGE fitsset.seed(4)ids <- IVBG %>%distinct(ID, .keep_all =TRUE) %>%group_by(area) %>%slice_sample(n =30)
Can we fit a single Sharpe Schoolfield with area-specific parameters with brms?
# Again, here are the data we are fitting:ggplot(dat, aes(temp, rate, color =factor(area, order$area))) +geom_point(data = dat, aes(temp, rate, color =factor(area, order$area)), size =0.6)
hist(dat$rate)
# Here's the equation (from the r package rTPC):# > sharpeschoolhigh_1981# function (temp, r_tref, e, eh, th, tref) # {# tref <- 273.15 + tref# k <- 8.62e-05# boltzmann.term <- r_tref * exp(e/k * (1/tref - 1/(temp + # 273.15)))# inactivation.term <- 1/(1 + exp(eh/k * (1/(th + 273.15) - # 1/(temp + 273.15))))# return(boltzmann.term * inactivation.term)# }# Add in fixed parametersdat$bk <-8.62e-05dat$tref <-8+273.15# (better visualization including bounds further down)n =10000hist(rnorm(0.8, 1, n = n), main ="e", xlim =c(0, 5))
hist(rnorm(2, 1, n = n), main ="eh", xlim =c(0, 7))
hist(rnorm(0.25, 0.5, n = n), main ="r_tref", xlim =c(0, 2.5))
hist(rnorm(12, 5, n = n), main ="th", xlim =c(5, 30))
# These work nicely with the pp_checkprior <-c(prior(normal(0.8, 1), nlpar ="e", lb =0),prior(normal(2, 1), nlpar ="eh", lb =0),prior(normal(0.3, 0.5), nlpar ="rtref", lb =0),prior(normal(10, 5), nlpar ="th", lb =0))# Now to a prior predictive checkfit_prior <-brm(bf(rate ~ (rtref *exp(e/bk * (1/tref -1/(temp +273.15)))) / (1+exp(eh/bk * (1/(th +273.15) -1/(temp +273.15)))), rtref + e + eh + th ~1, nl =TRUE),data = dat,cores =2,chains =2,iter =1500,sample_prior ="only",seed =9,prior = prior)
Compiling Stan program...
Trying to compile a simple C file
Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG -I"/Users/maxlindmark/Library/Caches/org.R-project.R/R/renv/cache/v5/R-4.2/x86_64-apple-darwin17.0/Rcpp/1.0.10/e749cae40fa9ef469b6050959517453c/Rcpp/include/" -I"/Users/maxlindmark/Dropbox/Max work/R/perch-growth/renv/library/R-4.2/x86_64-apple-darwin17.0/RcppEigen/include/" -I"/Users/maxlindmark/Dropbox/Max work/R/perch-growth/renv/library/R-4.2/x86_64-apple-darwin17.0/RcppEigen/include/unsupported" -I"/Users/maxlindmark/Dropbox/Max work/R/perch-growth/renv/library/R-4.2/x86_64-apple-darwin17.0/BH/include" -I"/Users/maxlindmark/Library/Caches/org.R-project.R/R/renv/cache/v5/R-4.2/x86_64-apple-darwin17.0/StanHeaders/2.26.28/35b5aee9ec9507aca2e021997a9e557e/StanHeaders/include/src/" -I"/Users/maxlindmark/Library/Caches/org.R-project.R/R/renv/cache/v5/R-4.2/x86_64-apple-darwin17.0/StanHeaders/2.26.28/35b5aee9ec9507aca2e021997a9e557e/StanHeaders/include/" -I"/Users/maxlindmark/Library/Caches/org.R-project.R/R/renv/cache/v5/R-4.2/x86_64-apple-darwin17.0/RcppParallel/5.1.7/a45594a00f5dbb073d5ec9f48592a08a/RcppParallel/include/" -I"/Users/maxlindmark/Library/Caches/org.R-project.R/R/renv/cache/v5/R-4.2/x86_64-apple-darwin17.0/rstan/2.32.3/4ac5b7639d28cd4fab19baaf46f33c6a/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DUSE_STANC3 -DSTRICT_R_HEADERS -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION -D_HAS_AUTO_PTR_ETC=0 -include '/Users/maxlindmark/Library/Caches/org.R-project.R/R/renv/cache/v5/R-4.2/x86_64-apple-darwin17.0/StanHeaders/2.26.28/35b5aee9ec9507aca2e021997a9e557e/StanHeaders/include/stan/math/prim/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -I/usr/local/include -fPIC -Wall -g -O2 -c foo.c -o foo.o
In file included from <built-in>:1:
In file included from /Users/maxlindmark/Library/Caches/org.R-project.R/R/renv/cache/v5/R-4.2/x86_64-apple-darwin17.0/StanHeaders/2.26.28/35b5aee9ec9507aca2e021997a9e557e/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
In file included from /Users/maxlindmark/Dropbox/Max work/R/perch-growth/renv/library/R-4.2/x86_64-apple-darwin17.0/RcppEigen/include/Eigen/Dense:1:
In file included from /Users/maxlindmark/Dropbox/Max work/R/perch-growth/renv/library/R-4.2/x86_64-apple-darwin17.0/RcppEigen/include/Eigen/Core:88:
/Users/maxlindmark/Dropbox/Max work/R/perch-growth/renv/library/R-4.2/x86_64-apple-darwin17.0/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
namespace Eigen {
^
/Users/maxlindmark/Dropbox/Max work/R/perch-growth/renv/library/R-4.2/x86_64-apple-darwin17.0/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
namespace Eigen {
^
;
In file included from <built-in>:1:
In file included from /Users/maxlindmark/Library/Caches/org.R-project.R/R/renv/cache/v5/R-4.2/x86_64-apple-darwin17.0/StanHeaders/2.26.28/35b5aee9ec9507aca2e021997a9e557e/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
In file included from /Users/maxlindmark/Dropbox/Max work/R/perch-growth/renv/library/R-4.2/x86_64-apple-darwin17.0/RcppEigen/include/Eigen/Dense:1:
/Users/maxlindmark/Dropbox/Max work/R/perch-growth/renv/library/R-4.2/x86_64-apple-darwin17.0/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
#include <complex>
^~~~~~~~~
3 errors generated.
make: *** [foo.o] Error 1
Start sampling
# Global prior predictive check in relation to data. Doesn't loo too informative...dat %>%data_grid(temp =seq_range(temp, n =51)) %>%ungroup() %>%mutate(bk =8.62e-05,tref =8+273.15) %>%add_epred_draws(fit_prior) %>%ggplot(aes(temp, y = .epred)) +stat_lineribbon(.width =c(0.75), alpha =0.3, color ="black", fill ="black") +geom_point(data = dat, aes(temp, rate, color =factor(area, order$area)), size =0.6) +scale_color_manual(values = colors, name ="Site") +labs(y ="Expectation of the posterior predictive distribution", x ="Temperature (°C)") +NULL
ungroup: no grouping variables
mutate: new variable 'bk' (double) with one unique value and 0% NA
new variable 'tref' (double) with one unique value and 0% NA
ggsave(paste0(home, "/figures/supp/prior_predictive_check.pdf" ), width =15, height =11, units ="cm")# Now make sure it converges with real data but no random effectsfit_global <-brm(bf(rate ~ (rtref *exp(e/bk * (1/tref -1/(temp +273.15)))) / (1+exp(eh/bk * (1/(th +273.15) -1/(temp +273.15)))), rtref + e + eh + th ~1,nl =TRUE),data = dat,cores =4,chains =4,iter =4000, seed =9,sample_prior ="yes",family = student,prior = prior )
Compiling Stan program...
Trying to compile a simple C file
Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/Rcpp/include/" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/unsupported" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/BH/include" -I"/Users/maxlindmark/Library/R/x86_64/4.2/library/StanHeaders/include/src/" -I"/Users/maxlindmark/Library/R/x86_64/4.2/library/StanHeaders/include/" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppParallel/include/" -I"/Users/maxlindmark/Library/R/x86_64/4.2/library/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DUSE_STANC3 -DSTRICT_R_HEADERS -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION -D_HAS_AUTO_PTR_ETC=0 -include '/Users/maxlindmark/Library/R/x86_64/4.2/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -I/usr/local/include -fPIC -Wall -g -O2 -c foo.c -o foo.o
In file included from <built-in>:1:
In file included from /Users/maxlindmark/Library/R/x86_64/4.2/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Dense:1:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Core:88:
/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
namespace Eigen {
^
/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
namespace Eigen {
^
;
In file included from <built-in>:1:
In file included from /Users/maxlindmark/Library/R/x86_64/4.2/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Dense:1:
/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
#include <complex>
^~~~~~~~~
3 errors generated.
make: *** [foo.o] Error 1
Start sampling
plot(fit_global)
pp_check(fit_global)
Using 10 posterior draws for ppc type 'dens_overlay' by default.
Ok, seems to work with priors and everything. They are roughly as broad as can be and still have a fit. Now fit the full model, with random area effects, more iterations and chains!
Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/Rcpp/include/" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/unsupported" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/BH/include" -I"/Users/maxlindmark/Library/R/x86_64/4.2/library/StanHeaders/include/src/" -I"/Users/maxlindmark/Library/R/x86_64/4.2/library/StanHeaders/include/" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppParallel/include/" -I"/Users/maxlindmark/Library/R/x86_64/4.2/library/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DUSE_STANC3 -DSTRICT_R_HEADERS -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION -D_HAS_AUTO_PTR_ETC=0 -include '/Users/maxlindmark/Library/R/x86_64/4.2/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -I/usr/local/include -fPIC -Wall -g -O2 -c foo.c -o foo.o
In file included from <built-in>:1:
In file included from /Users/maxlindmark/Library/R/x86_64/4.2/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Dense:1:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Core:88:
/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
namespace Eigen {
^
/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
namespace Eigen {
^
;
In file included from <built-in>:1:
In file included from /Users/maxlindmark/Library/R/x86_64/4.2/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Dense:1:
/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
#include <complex>
^~~~~~~~~
3 errors generated.
make: *** [foo.o] Error 1
Start sampling
print(fitbs, digits =4)
Family: student
Links: mu = identity; sigma = identity; nu = identity
Formula: rate ~ (rtref * exp(e/bk * (1/tref - 1/(temp + 273.15))))/(1 + exp(eh/bk * (1/(th + 273.15) - 1/(temp + 273.15))))
rtref ~ 1 + (1 | area)
e ~ 1 + (1 | area)
eh ~ 1 + (1 | area)
th ~ 1 + (1 | area)
Data: dat (Number of observations: 306)
Draws: 4 chains, each with iter = 5000; warmup = 2500; thin = 1;
total post-warmup draws = 10000
Group-Level Effects:
~area (Number of levels: 10)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
sd(rtref_Intercept) 0.0447 0.0384 0.0020 0.1434 1.0027 1649
sd(e_Intercept) 0.1823 0.1285 0.0089 0.4883 1.0009 2998
sd(eh_Intercept) 0.3051 0.3025 0.0083 1.0850 1.0031 1911
sd(th_Intercept) 0.8736 0.7670 0.0371 2.8345 1.0023 2302
Tail_ESS
sd(rtref_Intercept) 2800
sd(e_Intercept) 3025
sd(eh_Intercept) 4173
sd(th_Intercept) 4057
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
rtref_Intercept 0.4328 0.2079 0.2148 1.0089 1.0023 1418 3286
e_Intercept 0.8500 0.4512 0.1991 1.8886 1.0026 1442 3741
eh_Intercept 1.6155 0.4614 0.6903 2.5314 1.0004 3357 3018
th_Intercept 9.1672 3.8499 2.9851 17.2414 1.0021 1413 3481
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 0.0397 0.0020 0.0357 0.0435 1.0002 10449 7698
nu 24.6696 13.2991 8.2356 59.1964 1.0003 10512 7702
Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
Inspect the model
# Check fitlibrary(performance)
Attaching package: 'performance'
The following objects are masked from 'package:modelr':
mae, mse, rmse
fitbs
Family: student
Links: mu = identity; sigma = identity; nu = identity
Formula: rate ~ (rtref * exp(e/bk * (1/tref - 1/(temp + 273.15))))/(1 + exp(eh/bk * (1/(th + 273.15) - 1/(temp + 273.15))))
rtref ~ 1 + (1 | area)
e ~ 1 + (1 | area)
eh ~ 1 + (1 | area)
th ~ 1 + (1 | area)
Data: dat (Number of observations: 306)
Draws: 4 chains, each with iter = 5000; warmup = 2500; thin = 1;
total post-warmup draws = 10000
Group-Level Effects:
~area (Number of levels: 10)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(rtref_Intercept) 0.04 0.04 0.00 0.14 1.00 1649 2800
sd(e_Intercept) 0.18 0.13 0.01 0.49 1.00 2998 3025
sd(eh_Intercept) 0.31 0.30 0.01 1.08 1.00 1911 4173
sd(th_Intercept) 0.87 0.77 0.04 2.83 1.00 2302 4057
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
rtref_Intercept 0.43 0.21 0.21 1.01 1.00 1418 3286
e_Intercept 0.85 0.45 0.20 1.89 1.00 1442 3741
eh_Intercept 1.62 0.46 0.69 2.53 1.00 3357 3018
th_Intercept 9.17 3.85 2.99 17.24 1.00 1413 3481
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 0.04 0.00 0.04 0.04 1.00 10449 7698
nu 24.67 13.30 8.24 59.20 1.00 10512 7702
Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
r2_bayes(fitbs)
# Bayesian R2 with Compatibility Interval
Conditional R2: 0.373 (95% CI [0.301, 0.438])
r2_bayes(fit_global)
# Bayesian R2 with Compatibility Interval
Conditional R2: 0.120 (95% CI [0.057, 0.187])
Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/Rcpp/include/" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/unsupported" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/BH/include" -I"/Users/maxlindmark/Library/R/x86_64/4.2/library/StanHeaders/include/src/" -I"/Users/maxlindmark/Library/R/x86_64/4.2/library/StanHeaders/include/" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppParallel/include/" -I"/Users/maxlindmark/Library/R/x86_64/4.2/library/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DUSE_STANC3 -DSTRICT_R_HEADERS -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION -D_HAS_AUTO_PTR_ETC=0 -include '/Users/maxlindmark/Library/R/x86_64/4.2/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -I/usr/local/include -fPIC -Wall -g -O2 -c foo.c -o foo.o
In file included from <built-in>:1:
In file included from /Users/maxlindmark/Library/R/x86_64/4.2/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Dense:1:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Core:88:
/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
namespace Eigen {
^
/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
namespace Eigen {
^
;
In file included from <built-in>:1:
In file included from /Users/maxlindmark/Library/R/x86_64/4.2/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Dense:1:
/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
#include <complex>
^~~~~~~~~
3 errors generated.
make: *** [foo.o] Error 1
Start sampling
Warning: There were 1 divergent transitions after warmup. See
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.
Warning: Examine the pairs() plot to diagnose sampling problems
Warning: Dropping 'draws_df' class as required metadata was removed.
rename: renamed 4 variables (r_tref, e, eh, th)
pivot_longer: reorganized (r_tref, e, eh, th) into (parameter, value) [was 8000x4, now 32000x2]
mutate: new variable 'type' (character) with one unique value and 0% NA
Warning: Dropping 'draws_df' class as required metadata was removed.
rename: renamed 4 variables (r_tref, e, eh, th)
pivot_longer: reorganized (r_tref, e, eh, th) into (parameter, value) [was 8000x4, now 32000x2]
mutate: new variable 'type' (character) with one unique value and 0% NA
mutate: new variable 'iter' (integer) with one unique value and 0% NA
Compiling Stan program...
Trying to compile a simple C file
Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/Rcpp/include/" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/unsupported" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/BH/include" -I"/Users/maxlindmark/Library/R/x86_64/4.2/library/StanHeaders/include/src/" -I"/Users/maxlindmark/Library/R/x86_64/4.2/library/StanHeaders/include/" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppParallel/include/" -I"/Users/maxlindmark/Library/R/x86_64/4.2/library/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DUSE_STANC3 -DSTRICT_R_HEADERS -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION -D_HAS_AUTO_PTR_ETC=0 -include '/Users/maxlindmark/Library/R/x86_64/4.2/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -I/usr/local/include -fPIC -Wall -g -O2 -c foo.c -o foo.o
In file included from <built-in>:1:
In file included from /Users/maxlindmark/Library/R/x86_64/4.2/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Dense:1:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Core:88:
/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
namespace Eigen {
^
/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
namespace Eigen {
^
;
In file included from <built-in>:1:
In file included from /Users/maxlindmark/Library/R/x86_64/4.2/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Dense:1:
/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
#include <complex>
^~~~~~~~~
3 errors generated.
make: *** [foo.o] Error 1
Start sampling
Warning: There were 13 divergent transitions after warmup. See
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.
Warning: Examine the pairs() plot to diagnose sampling problems
Warning: Dropping 'draws_df' class as required metadata was removed.
rename: renamed 4 variables (r_tref, e, eh, th)
pivot_longer: reorganized (r_tref, e, eh, th) into (parameter, value) [was 8000x4, now 32000x2]
mutate: new variable 'type' (character) with one unique value and 0% NA
Warning: Dropping 'draws_df' class as required metadata was removed.
rename: renamed 4 variables (r_tref, e, eh, th)
pivot_longer: reorganized (r_tref, e, eh, th) into (parameter, value) [was 8000x4, now 32000x2]
mutate: new variable 'type' (character) with one unique value and 0% NA
mutate: new variable 'iter' (integer) with one unique value and 0% NA
Compiling Stan program...
Trying to compile a simple C file
Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/Rcpp/include/" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/unsupported" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/BH/include" -I"/Users/maxlindmark/Library/R/x86_64/4.2/library/StanHeaders/include/src/" -I"/Users/maxlindmark/Library/R/x86_64/4.2/library/StanHeaders/include/" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppParallel/include/" -I"/Users/maxlindmark/Library/R/x86_64/4.2/library/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DUSE_STANC3 -DSTRICT_R_HEADERS -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION -D_HAS_AUTO_PTR_ETC=0 -include '/Users/maxlindmark/Library/R/x86_64/4.2/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -I/usr/local/include -fPIC -Wall -g -O2 -c foo.c -o foo.o
In file included from <built-in>:1:
In file included from /Users/maxlindmark/Library/R/x86_64/4.2/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Dense:1:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Core:88:
/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
namespace Eigen {
^
/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
namespace Eigen {
^
;
In file included from <built-in>:1:
In file included from /Users/maxlindmark/Library/R/x86_64/4.2/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Dense:1:
/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
#include <complex>
^~~~~~~~~
3 errors generated.
make: *** [foo.o] Error 1
Start sampling
Warning: There were 4 divergent transitions after warmup. See
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.
Warning: Examine the pairs() plot to diagnose sampling problems
Warning: Dropping 'draws_df' class as required metadata was removed.
rename: renamed 4 variables (r_tref, e, eh, th)
pivot_longer: reorganized (r_tref, e, eh, th) into (parameter, value) [was 8000x4, now 32000x2]
mutate: new variable 'type' (character) with one unique value and 0% NA
Warning: Dropping 'draws_df' class as required metadata was removed.
rename: renamed 4 variables (r_tref, e, eh, th)
pivot_longer: reorganized (r_tref, e, eh, th) into (parameter, value) [was 8000x4, now 32000x2]
mutate: new variable 'type' (character) with one unique value and 0% NA
mutate: new variable 'iter' (integer) with one unique value and 0% NA
Compiling Stan program...
Trying to compile a simple C file
Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/Rcpp/include/" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/unsupported" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/BH/include" -I"/Users/maxlindmark/Library/R/x86_64/4.2/library/StanHeaders/include/src/" -I"/Users/maxlindmark/Library/R/x86_64/4.2/library/StanHeaders/include/" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppParallel/include/" -I"/Users/maxlindmark/Library/R/x86_64/4.2/library/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DUSE_STANC3 -DSTRICT_R_HEADERS -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION -D_HAS_AUTO_PTR_ETC=0 -include '/Users/maxlindmark/Library/R/x86_64/4.2/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -I/usr/local/include -fPIC -Wall -g -O2 -c foo.c -o foo.o
In file included from <built-in>:1:
In file included from /Users/maxlindmark/Library/R/x86_64/4.2/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Dense:1:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Core:88:
/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
namespace Eigen {
^
/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
namespace Eigen {
^
;
In file included from <built-in>:1:
In file included from /Users/maxlindmark/Library/R/x86_64/4.2/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Dense:1:
/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
#include <complex>
^~~~~~~~~
3 errors generated.
make: *** [foo.o] Error 1
Start sampling
Warning: Dropping 'draws_df' class as required metadata was removed.
rename: renamed 4 variables (r_tref, e, eh, th)
pivot_longer: reorganized (r_tref, e, eh, th) into (parameter, value) [was 8000x4, now 32000x2]
mutate: new variable 'type' (character) with one unique value and 0% NA
Warning: Dropping 'draws_df' class as required metadata was removed.
rename: renamed 4 variables (r_tref, e, eh, th)
pivot_longer: reorganized (r_tref, e, eh, th) into (parameter, value) [was 8000x4, now 32000x2]
mutate: new variable 'type' (character) with one unique value and 0% NA
mutate: new variable 'iter' (integer) with one unique value and 0% NA
Compiling Stan program...
Trying to compile a simple C file
Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/Rcpp/include/" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/unsupported" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/BH/include" -I"/Users/maxlindmark/Library/R/x86_64/4.2/library/StanHeaders/include/src/" -I"/Users/maxlindmark/Library/R/x86_64/4.2/library/StanHeaders/include/" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppParallel/include/" -I"/Users/maxlindmark/Library/R/x86_64/4.2/library/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DUSE_STANC3 -DSTRICT_R_HEADERS -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION -D_HAS_AUTO_PTR_ETC=0 -include '/Users/maxlindmark/Library/R/x86_64/4.2/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -I/usr/local/include -fPIC -Wall -g -O2 -c foo.c -o foo.o
In file included from <built-in>:1:
In file included from /Users/maxlindmark/Library/R/x86_64/4.2/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Dense:1:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Core:88:
/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
namespace Eigen {
^
/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
namespace Eigen {
^
;
In file included from <built-in>:1:
In file included from /Users/maxlindmark/Library/R/x86_64/4.2/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Dense:1:
/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
#include <complex>
^~~~~~~~~
3 errors generated.
make: *** [foo.o] Error 1
Start sampling
Warning: There were 1 divergent transitions after warmup. See
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.
Warning: Examine the pairs() plot to diagnose sampling problems
Warning: Dropping 'draws_df' class as required metadata was removed.
rename: renamed 4 variables (r_tref, e, eh, th)
pivot_longer: reorganized (r_tref, e, eh, th) into (parameter, value) [was 8000x4, now 32000x2]
mutate: new variable 'type' (character) with one unique value and 0% NA
Warning: Dropping 'draws_df' class as required metadata was removed.
rename: renamed 4 variables (r_tref, e, eh, th)
pivot_longer: reorganized (r_tref, e, eh, th) into (parameter, value) [was 8000x4, now 32000x2]
mutate: new variable 'type' (character) with one unique value and 0% NA
mutate: new variable 'iter' (integer) with one unique value and 0% NA
ts <- dat %>%data_grid(temp =seq_range(temp, n =100)) %>%mutate(bk =8.62e-05,tref =8+273.15) %>%add_epred_draws(fit_global, re_formula =NA, ndraws =8000)
mutate: new variable 'bk' (double) with one unique value and 0% NA
new variable 'tref' (double) with one unique value and 0% NA
p1 <-ggplot(ts) +geom_line(aes(temp, y = .epred, group = .draw), alpha = .1) +coord_cartesian(ylim =c(0.1, 0.4))p2 <-ggplot(ts) +geom_line(aes(temp, y = .epred, group = .draw), alpha = .1) +coord_cartesian(ylim =c(0.1, 0.4))p1 + p2
# Compute quantiles across spaghettiest_opt <- ts %>%group_by(.draw) %>%filter(.epred ==max(.epred)) %>%ungroup() %>%filter(!temp ==min(temp)) %>%# remove values where there was no optimumfilter(!temp ==max(temp))
# Make the main plot (conditional effect of temperature, with and without random effects)# Predictions without random effectsglob <- dat %>%data_grid(temp =seq_range(temp, n =100)) %>%mutate(bk =8.62e-05,tref =8+273.15) %>%add_epred_draws(fit_global, re_formula =NA) %>%#ungroup()
mutate: new variable 'bk' (double) with one unique value and 0% NA
new variable 'tref' (double) with one unique value and 0% NA
group_by: one grouping variable (area)
ungroup: no grouping variables
mutate: new variable 'bk' (double) with one unique value and 0% NA
new variable 'tref' (double) with one unique value and 0% NA
ungroup: no grouping variables
ggsave(paste0(home, "/figures/sharpe_school_bayes.pdf"), width =17, height =17, units ="cm")#ggsave(paste0(home, "/figures/for-talks/sharpe_school_bayes.pdf"), width = 14, height = 14, units = "cm")
Area-specific T_opts as a ridgeplot
# Calculate T_opt by areatsa <- dat %>%data_grid(area =unique(dat$area),temp =seq_range(temp, n =100)) %>%mutate(bk =8.62e-05,tref =8+273.15) %>%add_epred_draws(fitbs, re_formula =NULL, ndraws =10000)
mutate: new variable 'bk' (double) with one unique value and 0% NA
new variable 'tref' (double) with one unique value and 0% NA
t_opt_a <- tsa %>%group_by(.draw, area) %>%filter(.epred ==max(.epred)) %>%ungroup() %>%filter(!temp ==min(temp)) %>%# remove values where there was no optimumfilter(!temp ==max(temp))
mutate: new variable 'above' (character) with 2 unique values and 0% NA
group_by: one grouping variable (area)
count: now 16 rows and 3 columns, one group variable remaining (area)
pivot_wider: reorganized (above, n) into (N, Y) [was 16x3, now 10x3]
mutate (grouped): new variable 'prop' (double) with 7 unique values and 40% NA
# A tibble: 10 × 4
# Groups: area [10]
area N Y prop
<chr> <int> <int> <dbl>
1 FM 10 27 0.730
2 BT 7 15 0.682
3 SI_EK 14 10 0.417
4 JM 37 23 0.383
5 FB 44 3 0.0638
6 BS 16 1 0.0588
7 HO 29 NA NA
8 MU 18 NA NA
9 RA 14 NA NA
10 SI_HA NA 38 NA
# Without impacted?left_join(dat, t_opt_a_sum, by ="area") %>%filter(!area %in%c("BT", "SI_EK")) %>%filter(temp < median) %>%nrow()
left_join: added 3 columns (median, lwr, upr)
> rows only in x 0
> rows only in y ( 0)
> matched rows 306
> =====
> rows total 306
filter: removed 46 rows (15%), 260 rows remaining
filter: removed 92 rows (35%), 168 rows remaining
[1] 168
dat %>%group_by(area) %>%summarise(n =n()) %>%arrange(desc(n))
group_by: one grouping variable (area)
summarise: now 10 rows and 2 columns, ungrouped
# A tibble: 10 × 2
area n
<chr> <int>
1 JM 60
2 FB 47
3 SI_HA 38
4 FM 37
5 HO 29
6 SI_EK 24
7 BT 22
8 MU 18
9 BS 17
10 RA 14